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BENCHMARK COMPUTATIONS OF STRESSES IN A SPHERICAL 
DOME WITH SHELL FINITE ELEMENTS* 


ANTTI H. NIEMit 

Abstract. We present a computational framework for analysing thin shell structures using the 
finite element method. The framework is based on a mesh-dependent shell model which we derive 
from the general laws of three-dimensional elasticity. We apply the framework for the so called 
Girkmann benchmark problem involving a spherical shell stiffened with a foot ring. In particular, 
we compare the accuracy of different reduced strain four-node elements in this context. We conclude 
that the performance of the bilinear shell finite elements depends on the mesh quality but reasonable 
accuracy of the quantities of interest of the Girkmann problem can be attained in contrast to earlier 
results obtained with general shell elements for the problem. 
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1. Introduction. The role of thin shell theory in structural analysis has changed 
dramatically over centuries from center stage to supporting cast, partly because of 
the advent of the finite element method. This paper will bring shell theory back to 
the forefront by studying a family of four-node finite elements based on it. The study 
is carried out in the context of a challenging benchmark problem involving a stiffened 
doubly curved thin shell. 

Conventional shell finite element formulations involve various explicit and implicit 
modelling assumptions that extend beyond the limits of mathematical convergence 
theory currently available. The theoretical problems arise from the fact that the 
general shell elements used in industrial FEA have been developed through “finite 
element modelling” where the connection of the discretization to the actual differential 
operators of the mathematical shell model is obscure, cf. [30]. 

A prerequisite for traditional FE error analysis is a well-posed variational problem 
formulated in some suitable Hilbert space. For shells such mathematical models can 
be formulated using differential geometry of surfaces and theoretically stable formu¬ 
lations have been analyzed e.g. in [24, 1, 6]. However, these works deal only with 
the bending-dominated deformations and do not address the membrane-dominated 
and intermediate cases which are very important for instance in civil and structural 
engineering. 

On the other, it is possible to interpret the modelling assumptions of conventional 
shell elements in context of the mathematical models. For instance, the degenerated 
solid approach associated to certain four-node elements has been translated to explicit 
strain reduction procedure within a specific shell model in [17, 18] and that procedure 
has been numerically analysed in [11, 12, 22, 19]. This line of research is not limited to 
the theoretical analysis of existing formulations only. It can also be used to enhance 
the formulations and develop new ones as shown in [20], where a new four-node shell 
finite element of arbitrary quadrilateral shape was developed based on shell theory. 
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Thanks to modern computation technology, such as the /ip-adaptive finite ele¬ 
ment method [7, 33], shell analysis can also be based directly to three-dimensional 
elasticity theory. Such an approach rules out the modelling errors arising from the 
simplifications of dimensionally reduced structural models but requires more degrees 
of freedom for the discrete model. Also, if simplified representations of the stress 
state such as the stress resultants are needed, they must be post-processed from the 
three-dimensional stress field and this can be non-trivial. 

A model problem called the Girkmann problem, which was revived some time 
ago, highlights the above complications rather dramatically, see [31, 21, 28, 8]. The 
problem involves a concrete structure consisting of a spherical dome stiffened by a 
foot ring under a dead gravity load. The task is to determine the values of the shear 
force and the bending moment at the junction between the dome and the ring as well 
as the maximum bending moment in the dome. 

The problem was initially presented and solved analytically in the text book 
[10]. More recently, in the bulletin of the International Association of Computational 
Mechanics (lACM) [26], the problem was posed as a computational challenge to the 
finite element community. The purpose of the challenge was to find out how the 
process of verification, that is the process of building confidence that an approximative 
result is within a given tolerance of the exact solution to the mathematical model, 
is carried out by the community given the Girkmann problem. The results, that 
are summarized in [27, 31], without attribution and details on how verification was 
actually performed, are scandalous. Out of the 15 results submitted, 11 have a very 
large dispersion and are not within any acceptable tolerance of the reference values 
computed in [31, 21, 28] using different models and formulations. 

So far detailed verification studies have been published for the axisymmetric mod¬ 
els based on elasticity theory as well as axisymmetric dimensionally reduced models. 
In [31], the p-version of FEM was used in conjunction with the extraction procedure 
of [3] to compute accurate values for the quantities of interest. Similar approach with 
the hp-version of FEM was taken in [21], where also the axisymmetric h-version with 
selective reduced integration was successfully employed to discretize the dimensionally 
reduced model. 

In the present work, we introduce a finite element framework for thin shell analy¬ 
sis and use the Girkmann problem to demonstrate its possibilities. More precisely, we 
benchmark different variants of the MITC-type shell element proposed earlier by the 
author in [20] by modelling a quarter of the dome and by using symmetry boundary 
conditions. We show that the performance of the formulations varies depending on 
the performed strain reductions and the mesh regularity. Nevertheless, reasonable ac¬ 
curacy for the main quantities of interest is obtained contrary to the earlier published 
results obtained with general shell elements. 

The paper is organized as follows. In the next Section, we develop the shell theory 
used to construct the different FE methods. The finite element methods are described 
in Section 3 together with a discussion on their well-posedness and implementation 
aspects. Section 4 is devoted to description of the analysis procedure for the Girkmann 
problem and numerical results. The paper ends with conclusions in Section 5. 

2. Variational Formulation of the Shell Problem. We will employ the Ein¬ 
stein summation convention so that Greek indices range over the values 1, 2 while 
Latin indices have the values 1, 2, 3. The former refer to surface coordinates while 
the latter refer to general three-dimensional curvilinear coordinates. Our surface 
coordinate systems can be assumed orthonormal so that we can formulate our strain- 
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displacement relations and constitutive laws directly in terms of physical components 
and traditional partial derivatives which will be denoted by a comma. Moreover, 
Euclidean vectors are displayed using overhead arrows, whereas boldface notation 
is reserved for the column vectors and matrices storing the components of different 
surface tensors in the assumed orthonormal coordinate system. 

A shell domain C of constant thickness t is defined as 

(2.1) n= ${K X i-t/2,t/2)), 
where the shell mapping (p is of the form 

(2.2) ${x,yX)=rix,y) + Cnix,y). 

Here the parametric surface r{x, y) represents the middle surface of the shell and 
n{x, y) is the unit normal vector to the middle surface. Thus = cc, ^2 = y, and 
^3 = C constitute a curvilinear coordinate system in 3-space, called shell coordinates. 
In the following, we imagine that the middle surface is defined as 

(2.3) f{x,y) = xii+yi 2 +fix,y)i 3 , {x,y) G K, 

where ii, ^ 2 , *3 are fixed Cartesian unit vectors. Moreover, we assume that the middle 
surface differs only little from the coordinate plane K, i.e. the shell is shallow. More 
precisely, we assume that the function / : AT —>■ R is smooth in the curvature length 
scale R defined by 

K oi,p 

If hjc = diam(Ar) and R is taken as the length unit, the shallowness assumption can 
be formulated as 

(2.4) /,„ = 0{hK), 
where h^ = h^/R < I. 

Under the shallowness assumption (2.4), the tangent basis vectors 

(2.5) ea{x,y) = ia + f,aix,y)i3 

associated to the middle surface parametrization (2.3), are orthonormal within the 
accuracy of 0{h\). 

2.1. Shell Kinematics. According to the standard kinematic hypothesis we 
assume that the displacement vector can be written in the form 

(2.6) U{x,y,C) = {ux{x,y) + C9x{x,y))ex{x,y) + w{x,y)n{x,y), 

where u = {ui,U 2 ) are the tangential displacements of the middle surface, w is the 
transverse deflection, and the quantities 9 = { 61 , 62 ) are the angles of rotation of the 
normal. 

Referring to the curvilinear coordinates = x, ^2 = y, and ^3 = C,, the linearized 
Green-Lagrange strain tensor is defined by 

Gi = ^ ■ U,i): i, j = 1, 2,3. 


(2.7) 
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We find directly from (2.6) that 

(2.8) = {u\,a + 6^A,a)eA + (ma + i?.d\)ex,a + W,an + Wn a 
for a = 1, 2 and 

(2.9) f7,3 = ^ASA. 

Similarly, combination of (2.2), (2.3) and (2.5) yields 

(2.10) = Ca + Cfl.a, 0=1,2, 
and 

( 2 . 11 ) ^,3 = n - 

The components of the strain tensor are represented as power series of the variable 
C,. If we take into account two terms, the in-plane strains may be written as 

(2.12) Ca/3 ~ T 

Using relations (2.10) and (2.8), the membrane strain tensor Eq,^, which arises from 
stretching of the deformed middle surface, can be written as 

(2.13) £^a/3 ~ ^a/3^, 

where 

(2.14) ba/3 = -e; ■ n,p = , ^ , a,/3 = l,2, 

\/i + /a/a 

are the coefficients of the second fundamental form of the middle surface. In (2.13), 
terms multiplied by the ‘rotation coefficients’ Ca ■ e\^p have been neglected as quanti¬ 
ties of relative order 0{hK) based on the shallowness assumption (2.4). We content 
ourselves here to first order of accuracy since, in general, the coefficients cannot 
be approximated more precisely with linear or bilinear interpolation functions, see 
Section 3.2. 

Introducing the coefficients of the third fundamental form of the middle surface 
Ca/3 — ‘ ^,/3, tT, — 1,2, 

the elastic curvature tensor Kap, which arises from bending of the deformed middle 
surface, comes out as 

(2.15) ^q/ 3 ~ ^(^a,/3 T ^/3,a) '^{bo^xUx p ^/3A^A,a), CT, /5 — 1,2 

based again directly on (2.10) and (2.8). Here, terms multiplied by Ca-ex^p or fia-ex^p 
have been neglected as quantities of relative order 0 {hK)- 

It is possible to simplify the bending strain expressions by sacrificing their tenso- 
rial invariance. It is straightforward to verify that 
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within the adopted accuracy, so that we may write (2.15) component-wise as 


Kii Ri 6*1,1 + hi2{hi2W - U2,i ) - fciieii, 

K22 ~ ^2,2 + hi2{bi2W — ^1,2) — 6'22£22, 

«^12 ~ 2 + ^2,1 + bii{hi2W — ^1,2) -I- b22{hl2W — M2,l))-^22)- 


In these expressions, the contribution of the terms buEu, 622222 and 612(211 -I- 222 ) to 
the maximum in-plane strains at the outer and inner surfaces of the shell is of relative 
order 0{t/R) only. Therefore, the number of terms in the kinematic relations can be 
slightly reduced by retaining only the underlined terms in the calculations. 

Finally, the transverse shear strains are defined as 


(2.16) 




These can be written in terms of the displacements by first noting that since Ca and 
n are orthogonal, we have ea,j3 ■ n = —Ca ■ n,p, and consequently bap = n ■ ea,p- Now, 
combination of (2.11) with (2.8) and (2.10) with (2.9) according to (2.7) yields 

(2.17) 'ya=0a + b\aU\+W^a, 0=1,2, 


and completes the description of shell kinematics. 

2.2. Constitutive Equations. Assuming linearly elastic isotropic material with 
Poisson ratio v and Young modulus A, the constitutive law relating stresses and 
strains can be written with respect to the approximately orthogonal shell coordinate 
system (a:, y, Q attached to the middle surface as 


^a/3 — . 2 [(^ ^)2a/3 T ^^XX^aj3\ ; 

(2.18) 

Tck3 ~ \ 2ck3, O, /3 = 1, 2, 

l + V 

where dap is the Kronecker delta. The elastic coefficients in the above formula have 
been modified to yield so called plane stress state tangent to the middle surface. 
This modification is necessary to avoid Poisson locking in context of the kinematic 
assumption ( 2 . 6 ). 

We follow the standard convention of structural mechanics and introduce the 
internal forces and moments which are the stress resultants and stress couples per 
unit length of the middle surface. These can now be defined as 



rt/2 

rt/2 

rt/2 

(2.19) 

'^a/3 — / ^a/3 — 

1 Qoi — 

/ CTaS dC, 


J~t/2 

f-t/2 

l-t/2 


and correspond to the membrane forces, bending moments and transverse shear forces 
in static equilibrium considerations. 

2.3. Potential Energy Functional. Integrals over the shell domain can be 
evaluated in terms of the assumed shell coordinates as 

f {-IdriR::! f [ {■}dC,dxdy. 

Jn Jk J-t /2 
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Combination of (2.19), (2.18), (2.16) and (2.12) allows us to write the elastic strain 
energy functional as 


( 2 . 20 ) {7/f(M,W,0) = i [ a^lSea/B dfl Ki I- f {UapSap + qala + mapKap) dxdy 

^ JQ, JK 


where 


( 2 . 21 ) 


Et 


a/3 — 2 ^)^a/3 “1” j 


Qa = 


1 — E 

Et 


:1a 


2(1 + ix) 

Et^ 

ma/3 = J^2(l - 1,2^ “ ’^)Kal3 + l^KxxSap] 


and the strains are given in terms of the displacements in (2.13), (2.17) and (2.15). 

Similarly, the potential energy corresponding to external distributed surface forces 
(/i,/ 2 ,p) and moments (ri,r 2 ) is 


( 2 . 22 ) 


Vk{u,w, 



{fxux +pw + TxOx) dxdy 


and the total energy is given by the sum 


(2.23) Ek{u, w, 0) = Uk{u, w, 6 ) + Vk{u, w, 6 ). 

3. Finite Element Methods. We assume that the whole shell domain fl is 
formed as a union of domains of the form (2.1) as 

^ = U 

KeCh 

where Ch stands for a mesh of convex quadrilaterals K corresponding to parametriza- 
tions of patches of the shell middle surface according to (2.3). We also assume that 
the whole middle surface 


5= U tk{K) 

KcCh 

is a smooth surface and that it can be described alternatively by a single, global 
parametrization It follows that the transformations between the local and 

global coordinate systems Tk = op, K G Ch, are diffeomorphisms. 

Without losing generality, we may assume that the coordinates ^ 1,^2 are isother¬ 
mal. If {5i(Ci)^ 2 ) 552 ('Ci 5 C 2 )} are the corresponding orthonormal tangent vectors of 
the middle surface, and Ua, w, and 9a stand for the associated displacement and 
rotation components, then these components are related to the local components Ua , 
and 9^ as 

(3.1) Ua oTk = Uxgx ■ la, oTk = W, 9^ oTk= 9xgx ■ la 

according to (2.5) and (2.6). The total potential energy of the structure is expressed 
as the sum of element-wise contributions such that 

E{u,w,e)^ EK{u^,w^,e^), 

K&Ch 


(3.2) 
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where {u,w,6) is the globally defined generalized displacement field. The solution 
of the problem is determined according to the principle of minimum potential energy 
from the condition 


E{u,w,6)= min E{v,z,^), 


where the energy space U is defined as the set of those kinematically admissible 
generalized displacement fields for which the energy functional is finite. The existence 
of a unique minimizer (for max^f sufficiently small) follows from the well-posedness 
of the corresponding Reissner-Naghdi shell model, see e.g. [5]. 

It is now straightforward to formulate a finite element method where each dis¬ 
placement component is approximated separately as in the space 

(3.3) Uh = {iu,w,e)GU : (m^, , 6^) e [Qi{K)f VK G Ch}, 

where Qi{K) denotes the standard space of isoparametric bilinear functions on K 
and ,9^) is the local generalized displacement field defined in (3.1). 

3.1. Strain Reduction Techniques. To avoid locking when approximating 
bending-dominated problems, membrane and transverse shear strains must be re¬ 
duced. To introduce the different methods, we denote by Fk = (xk, Uk) the bilinear 
mapping of the reference square K = [—1,1] x [—1,1] onto K and by 

( dxK dxK \ 

Si Si] 

dx dy J 

the Jacobian matrix of Fpc- Here {x,y) are the coordinates on K. 

We start by defining on the reference square K the function spaces 

(3.4) S{K) = {s=(^^+f^ a,b,c,dGR} 

and 

(3.5) MiK) = {r=(^^+J^y : a,b,c,d,eGR} 

for the reduced transverse shear strains and membrane strains, respectively. The 
canonical degrees of freedom associated with S{K) are 


(3.6) s i-A / s^fds for every edge e of RT, 

J e 

whereas the degrees of freedom associated with M{K) are defined as 


zT 


(3.7) 


r i-A i rids for every edge e of K, 


r I—s- / ti 2 dxdy. 

J k 
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The corresponding spaces associated to a general element K € Ch are then defined 
using covariant transformation formulas as 

(3.8) S{K) = {s = {J-/s) o = Sk{s) : a € S{k)} 
and 

(3.9) M{K) = {t = r/0 o = MK(r) : f G M{k)} 

In case of general quadrilateral elements with non-constant Jacobian matrices, the 
transformation in (3.9) must be fixed to a single orientation e.g. at the midpoint of 
K so that J = J(0, 0). Otherwise the function space (3.9) for the reduced membrane 
strains and stresses does not necessarily include constant fields which would degrade 
the accuracy of the formulation. 

Denoting hy A: H^{k) —>■ S{k) and 11 ^ : H^{k) -G M{k) the interpolation 
operators associated to the degrees of freedom (3.6) and (3.7), the corresponding 
projectors for a general K G Ch are defined as 

TAk = -Mk o o and Ax = Sx o o 5^^ 

The transformation rule (3.8) guarantees that the degrees of freedom (3.6) are pre¬ 
served on K\ 


si-G s^t ds = 0 for every edge e of K, 

J e 

but the analogous statement does not hold in context of (3.9) for general meshes. 

The finite element introduced in (3.4) and (3.6) is a well known edge element 
denoted by the symbol RTcf in the recently introduced Periodic Table of the Finite 
Elements [2] and has been described in the context of plate bending e.g. in [13, 4]. 
The element (3.4), (3.7) is less customary at least in the mathematical literature - 
probably because of its non-elegant extension to quadrilateral shapes. In any case 
the historical roots of the element are very deep in the literature on finite element 
technology for plane elasticity. Our current formulation corresponds essentially to 
the stress field of the Pian-Sumihara element introduced in [23] which in turn may 
be viewed as an extension of the nonconforming displacement methods introduced by 
Wilson et. al. [35] and Turner et. al. [34]. We refer the reader to [25] for the complete 
mathematical theory of these formulations in context of plane elasticity. In the present 
context, the convergence theory is confined to special cases involving cylindrical or 
globally shallow shells on special meshes, see [11, 12, 22, 19]. 

We shall use the label MITC4C for the formulation for which only the transverse 
shear strains are projected into the space (3.8) 

(3.10) 7 ^ Axj 

and the label MITC4S for the formulation where also the membrane strains are pro¬ 
jected: 

(3.11) e ^ Il/i-e, 7 ^ Ax'j 

when evaluating the strain energy according to (2.20) and (2.21). 
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In addition, we consider stabilized variants of both methods, where the shear 
modulus G = 2 {^„) (2-21) is modified as 

(3.12) G -A Gk = ,2 ,2 G. 

Here aK is a positive stabilization parameter independent of t and Hk- This stabiliza¬ 
tion idea originates from the corresponding plate bending elements, see e.g. [32, 15]. 

Finally, the abbreviation D1SP4 is used for the standard displacement method 
without any strain reduction or stabilization. 

3.2. Implementation. The implementation of the present formulation follows 
the standard steps used to construct quadrilateral plane-elastic elements with some 
twists. The shape functions are the usual bilinear ones on the reference square K and 
the computation of derivatives and numerical integration using the Gauss quadrature 
are performed in the canonical way. 

The hrst twist concerns determination of the physical element K from the general 
surface mesh consisting of quadrilateral elements. As it may happen, that the four 
nodes of an element do not lie in the same plane, a straightening operation is needed 
in order to define the plane element K. This can be accomplished in many ways, but 
we follow the procedure of [16] described also in [20] that yields naturally also the 
directions of the local coordinate axes ii and * 2 . 

The second difference is related the strain-displacement relations (2.13),(2.15) and 
(2.17), which involve the coefficients of the second fundamental form in addition to 
the standard shape function derivatives. Within the limits of the local shallowness 
assumption, these coefficients can computed using the interpolated normal vector Uh 
as 


^a/3 ~ ioL * Ct;, (3 - 1; 2. 

To carry out the interpolation, the nodal normals are needed as geometric input data 
in addition to the node positions. The former are also used to construct the two 
orthonormal tangent direction gi , 52 used to enforce the continuity of the tangential 
displacements and normal rotations according to (3.1). 

The current implementation (assembly, solution and post-processing) is carried 
out using Mathematica while Gmsh is used for the mesh generation [14, 9]. An 
additional pre-processing step is the determination of the nodal normals that can be 
performed by using the analytic surface representation if available, or by averaging 
the normals of the elements sharing a common node. 

4. Numerical Results. We start by recalling the statement of the Girkmann 
benchmark problem from [26, 31]. The problem involves a concrete structure consist¬ 
ing of a spherical dome stiffened by a foot ring under a dead gravity load, see Fig. 1. 
The task is to determine the values of the transverse shear force and the meridional 
bending moment at the junction between the dome and the ring as well as the max¬ 
imum value of the bending moment in the dome assuming that the gravity load is 
equilibrated by a uniform pressure acting at the base of the ring. The material of the 
concrete is assumed to be linearly elastic, homogeneous and isotropic with vanishing 
Poisson ratio. The value of the Young modulus is specified as E = 20.59 GPa although 
it has no effect on the values of the quantities of interest. 

We follow here the classical approach, where the unknown reactions are taken to 
be the horizontal force R and the bending moment M and are assumed to be positive 
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Fig. 1. The Girkmann problem. Cross-section of the structure. 


when acting on the shell. In this splitting (shown in Fig. 1), the normal force N 
becomes determined from the vertical force balance as 


(4.1) 


N = 


-grp 
1 + cos a ’ 


where g = Ft is the vertical surface load density corresponding to the assumed weight 
density F = 32690 N/m^, a is the opening angle of the dome and tq its radius (Fig. 1). 
The shear force requested in the problem statement is then defined as Q = R/ sin a. 

The classical solution procedure can then be formulated as follows. Let us assume 
that A and T denote the horizontal displacement of the midpoint of the junction and 
the angle of rotation of the junction line, respectively. Then, the linearity of the 
material model implies that the following relations must hold 

EA = EA^ + kf^R + kf^M = EA^ + kfj^R + kf2M, 

E'S = FIT® + + k^^M = E'S^ + kfj^R + kf^M, 


where and denote the displacement and the rotation of the Shell/Ring due 
to known loads only and hj = 1,2 are the compliance constants associated to 
the unknown reactions R and M. These parameters can be defined separately for the 
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Fig. 2. Regularly refined mesh sequence. 





Fig. 3. Mesh sequence generated with the frontal method of Gmsh. 


shell and the ring by analysing sequentially the following load cases 

(Case 1) Gravity load (shell), equilibrating pressure (ring) and N as in (4.1), 
(Case 2) i? = lN/m, 

(Case 3) M = INm/m, 

and recording the values of the horizontal displacement and the rotation in each case. 

4.1. Convergence studies. Our focus is on the performance of the shell ele¬ 
ments so that we consider first the convergence of the parameters in (4.2) for the 
dome using two different kinds of mesh sequences with the maximum element size h 
approaching zero in both cases. The first mesh sequence is based on regular refine¬ 
ment of the initial mesh with three elements as shown in Fig. 2. The second mesh 
sequence is shown in Fig. 3 and is generated using the frontal algorithm of Gmsh 
restricted so that each edge has a fixed number of elements [29]. 

The results are shown in Tables 1 and 2 for the different mesh sequences. The 
reported values are normalized against the reference values 

EA^ ss -2.300E6N/m, fcfi ss 8.345 e 3, fcfa - 1.477E41/m, 

ss-9.338E5N/m2, fcl) ss-1.477 e 4 1/m, fcf 2 ~-5-113e 4 1/m^. 

computed using the same axisymmetric shell model as in [21]. The values are nearly 
identical to those of [28] computed using a classical shell model which neglects trans- 























Fig. 4. Total displacement in meters in the membrane-dominated load Case 1. MITC 4 S (left) 
versus MITCfC (right) on the frontal mesh with N = 32. 


verse shear deformations.^ 

The discretization parameter N stands for the number of elements per edge and 
is varied from 8 to 256. The required displacement and rotation values are calculated 
as the averages of the corresponding nodal values over the junction line and, for the 
stabilized variants of the methods, the value a = 0.2 is used for the stabilization pa¬ 
rameter. These results do not exhibit significant differences between the formulations 
apart from the standard displacement method DISP4 which locks in Cases 2 and 3 as 
expected. 

However, a more detailed investigation reveals a rather drastic difference between 
the different formulations in the solutions of the membrane-dominated Case 1. Fig. 4 
shows the total displacement of the shell mid-surface calculated on the frontal mesh 
with 16 elements per edge using the MITC4S and MITC4C formulations. It is evident 
that the MITC4S formulation suffers from numerical instabilities associated to the 
consistency error arising from the reduction of the membrane strains, cf. [22]. The 
instability is manifested here by the loss of symmetry of the numerical solution. It is 
intriguing that the circumferentially averaged displacements still converge and that 
the phenomenon disappears when a regular mesh is used as shown in Fig. 5. 

Finally, it should be pointed out that similar instabilities as shown in Fig. 4 
(left) are featured also by the lowest-order linear or bilinear shell elements employed 
currently in many industrial FEA programs. 

4.2. Section force and moment computations. In order to determine the 
unknown reactions R and M and the deformation of the stiffened shell, we need the 
compliance coefficients of the ring. These can be approximated by using the principle 
of virtual work as 

EK^ « 1.363E7N/m, « -2683, ~ 8418 1/m, 

-6.949E6N/m^ -84181/m, ~ 3.696 e 4 l/m^. 

These values are obtained by taking the aforementioned displacements A and as the 
only degrees of freedom as in [28]. This corresponds to the kinematic assumption that 
the ring cross section deforms as a rigid body. The associated 2x2 stiffness matrix 


^The unit of force adopted in reference [28] is Girkmann’s kilopond, 1 G = 9.807 N. 
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EA^o 


N 

DISP4 

MITC4C 

MITC4S 

Stab. MITC4C 

Stab. MITC4C 

8 

0.98 

0.70 

0.75 

0.68 

0.75 

16 

1.00 

0.89 

0.90 

0.89 

0.90 

32 

1.00 

0.97 

0.97 

0.97 

0.97 

64 

1.00 

0.99 

0.99 

0.99 

0.99 

128 

1.00 

1.00 

1.00 

1.00 

1.00 

256 

1.00 

1.00 

1.00 

1.00 

1.00 







N 

DISP4 

MITC4C 

MITC4S 

Stab. MITC4C 

Stab. MITC4S 

8 

1.00 

1.96 

2.19 

1.93 

2.12 

16 

0.99 

1.42 

1.47 

1.41 

1.47 

32 

1.00 

1.13 

1.14 

1.13 

1.14 

64 

1.01 

1.03 

1.04 

1.03 

1.04 

128 

1.01 

1.01 

1.01 

1.01 

1.01 

256 

1.01 

1.00 

1.00 

1.00 

1.00 




uS 

^11 



N 

DISP4 

MITC4C 

MITC4S 

Stab. MITC4C 

Stab. MITC4S 

8 

0.20 

0.51 

0.64 

0.52 

0.66 

16 

0.28 

0.78 

0.87 

0.78 

0.89 

32 

0.40 

0.93 

0.96 

0.93 

0.97 

64 

0.56 

0.98 

0.99 

0.98 

0.99 

128 

0.74 

1.00 

1.00 

1.00 

1.00 

256 

0.89 

1.00 

1.00 

1.00 

1.00 




uS 

^12 



N 

DISP4 

MITC4C 

MITC4S 

Stab. MITC4C 

Stab. MITC4S 

8 

0.04 

0.39 

0.56 

0.36 

0.53 

16 

0.09 

0.71 

0.83 

0.70 

0.83 

32 

0.17 

0.91 

0.95 

0.90 

0.95 

64 

0.32 

0.97 

0.99 

0.97 

0.99 

128 

0.56 

0.99 

1.00 

0.99 

1.00 

256 

0.80 

1.00 

1.00 

1.00 

1.00 




uS 

^22 



N 

DISP4 

MITC4C 

MITC4S 

Stab. MITC4C 

Stab. MITC4S 

8 

0.01 

0.52 

0.66 

0.67 

0.79 

16 

0.03 

0.78 

0.88 

0.81 

0.90 

32 

0.07 

0.93 

0.97 

0.94 

0.97 

64 

0.18 

0.98 

0.99 

0.98 

0.99 

128 

0.42 

1.00 

1.00 

1.00 

1.00 

256 

0.72 

1.00 

1.00 

1.00 

1.00 


Table 1 

Convergence of the compliance coefficients for the dome with respect to the mesh parameter N 
on the regular mesh sequence (Fig. 2). 
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EAl 


N 

DISP4 

MITC4C 

MITC4S 

Stab. MITC4C 

Stab. MITC4C 

8 

1.00 

0.74 

0.92 

0.72 

0.98 

16 

1.02 

0.91 

0.94 

0.90 

0.94 

32 

1.01 

0.97 

0.98 

0.97 

0.98 

64 

1.01 

0.99 

1.00 

0.99 

1.00 

128 

1.00 

1.00 

1.00 

1.00 

1.00 

256 

1.00 

1.00 

1.00 

1.00 

1.00 







N 

DISP4 

MITC4C 

MITC4S 

Stab. MITC4C 

Stab. MITC4S 

8 

0.96 

1.78 

1.83 

1.83 

1.81 

16 

0.96 

1.37 

1.45 

1.37 

1.45 

32 

0.98 

1.12 

1.05 

1.12 

1.05 

64 

0.99 

1.03 

1.01 

1.03 

1.01 

128 

1.00 

1.01 

1.00 

1.01 

1.00 

256 

1.00 

1.00 

1.00 

1.00 

1.00 




uS 

^11 



N 

DISP4 

MITC4C 

MITC4S 

Stab. MITC4C 

Stab. MITC4S 

8 

0.20 

0.48 

0.57 

0.50 

0.68 

16 

0.28 

0.71 

0.78 

0.75 

0.82 

32 

0.40 

0.89 

0.91 

0.91 

0.94 

64 

0.54 

0.96 

0.98 

0.97 

0.98 

128 

0.73 

0.99 

0.99 

0.99 

1.00 

256 

0.86 

1.00 

1.00 

1.00 

1.00 




uS 

^12 



N 

DISP4 

MITC4C 

MITC4S 

Stab. MITC4C 

Stab. MITC4S 

8 

0.04 

0.32 

0.42 

0.33 

0.46 

16 

0.08 

0.62 

0.70 

0.63 

0.72 

32 

0.17 

0.85 

0.88 

0.87 

0.90 

64 

0.30 

0.95 

0.97 

0.96 

0.97 

128 

0.53 

0.99 

0.99 

0.99 

0.99 

256 

0.75 

1.00 

1.00 

1.00 

1.00 




uS 

^22 



N 

DISP4 

MITC4C 

MITC4S 

Stab. MITC4C 

Stab. MITC4S 

8 

0.01 

0.34 

0.41 

0.69 

0.78 

16 

0.02 

0.70 

0.76 

0.77 

0.84 

32 

0.06 

0.88 

0.90 

0.91 

0.94 

64 

0.16 

0.96 

0.97 

0.97 

0.98 

128 

0.37 

0.99 

0.99 

0.99 

0.99 

256 

0.64 

1.00 

1.00 

1.00 

1.00 


Table 2 

Convergence of the compliance coefficients for the dome with respect to the mesh parameter N 
on the frontal mesh sequence (Fig. 3). 













Fig. 5. Total displacement in meters in the membrane-dominated load Case 1. MITC 4 S (left) 
versus MITCfC (right) on the regular mesh with N = 32. 


is evaluated by using numerical integration up to the machine precision in cylindri¬ 
cal coordinates over the exact pentagonal shape of the ring. The ring is assumed 
weightless here as in the original treatment by Girkmann and in the contemporary 
verification challenge. 

Substitution of (4.4) and (4.3) into (4.2) yields the values of the unknown hori¬ 
zontal section force and bending moment: 

(4.5) i?Ril467N/m & M «-37.36 Nm/m. 

Finally, the deformation of the stiffened dome can be computed by solving the shell 
problem with the active loads as in 

(Case 4) Gravity load, N according to (4.1) and R,M according to (4.5). 

In order to demonstrate the influence of the stabilization technique (3.12), we 
show in Figs. 6-7 the distribution of the meridional bending moment as computed 
with the different formulations along the left and right edges of the computational 
domain. The post-processing is carried out directly from the nodal rotations and only 
the range (20°, 40°) is shown here since the bending effects are confined to a narrow 
region near the edge. 

In this case, there is no visible difference between the MITC4S and MITC4C 
formulations. For instance, on the frontal mesh with N = 32, both formulations 
feature non-physical oscillations but the stabilization technique (3.12) improves the 
results as shown in Figs. 6 and 7. On the other hand, a feasible solution is obtained 
even without the stabilization when the mesh is regular as shown in Fig. 8. 

5. Conclusions. We have presented a finite element framework for analysing 
the structural response of thin elastic shells. The element stiffness matrices are com¬ 
puted according to shell theory taking into account transverse shear deformations. 
The framework enables enables explicit reduction of the membrane strains and the 
transverse shear strains in order to resolve numerical locking problems. 

Different variants of quadrilateral elements have been examined in the Girkmann 
problem and a detailed verification study has been performed. We have demonstrated 
that the MITC4S element with reduced membrane strains may feature numerical 
instabilities in membrane-dominated situations if the mesh is irregular. The MITG4G 
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Bending Moment 
Nm/m 



Bending Moment 
Nm/m 



Fig. 6. Distribution of the meridional bending moment in the stiffened dome calculated along 
the left and right edges of the computational domain, respectively. MITC 4 S formulation on the 
frontal mesh with N = 32. 


Bending Moment 
Nm/m 



Bending Moment 
Nm/m 



Fig. 7. Distribution of the meridional bending moment in the stiffened dome calculated along 
the left and right edges of the computational domain, respectively. Stabilized MITC 4 S formulation 
on the frontal mesh with N = 32. 


where only the transverse shear strains are reduced, does not have this drawback. 
On the other hand, the bending moment in the Girkmann problem can be calculated 
accurately with both methods on regular meshes. On irregular meshes some non¬ 
physical oscillations occur but these can be avoided by employing the transverse shear 
balancing as a stabilization technique. 
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